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Abstract 

By means of a three-dimensional amphiphilic lattice-Boltzmann model with short-range interactions for the description 
of ternary amphiphilic fluids, we study how the phase separation kinetics of a symmetric binary immiscible fluid is 
altered by the presence of the amphiphilic species. We find that a gradual increase in amphiphile concentration slows 
down domain growth, initially from algebraic, to logarithmic temporal dependence, and, at higher concentrations, from 
logarithmic to stretched-exponential form. In growth-arrested stretched-exponential regimes, at late times we observe the 
self-assembly of sponge mesophases and gyroid liquid crystalline cubic mesophases, hence confirming that (a) amphiphile- 
amphiphile interactions need not be long-ranged in order for periodically modulated structures to arise in a dynamics 
of competing interactions, and (b) a chemically-specific model of the amphiphile is not required for the self-assembly of 
cubic mesophases, contradicting claims in the literature. We also observe a structural order-disorder transition between 
sponge and gyroid phases driven by amphiphile concentration alone or, independently, by the amphiphile-amphiphile 
and the amphiphile-binary fluid coupling parameters. For the growth-arrested mesophases, we also observe temporal 
oscillations in the structure function at all length scales; most of the wavenumbers show slow decay, and long-term 
stationarity or growth for the others. We ascribe this behaviour to a combination of complex amphiphile dynamics 
leading to Marangoni flows. 

PACS numbers: mm.nn.xx 



* n.gonzalez-segredoOucl . ac .uk. Also at Departament de Fi'sica, Universitat Autonoma de Barcelona, 08193 

Bellaterra, Barcelona, Spain, 
t p . V . coveneyOucl .ac.uk 



1 



I. INTRODUCTION 



The term amphiphilic fluid is broadly used to denote multiphase fluids in which at least one species is of 
a surfactant nature. A surfactant molecule (from surface active agent, which we shall also refer to as an 
amphiphile) contains a polar headgroup attached to a hydrocarbon tail which, dispersed in a binary immiscible 
fluid mixture, such as oil and water, is driven towards and adsorbed at the interface between the two fluids. 
The selective chemical affinity between each part of the surfactant molecule and the components of the binary 
mixture is the mechanism responsible for such a taxis (jyj . Not only are amphiphilic fluids important in physical 
chemistry, structural biology, soft matter physics and materials science from a fundamental perspective, but 
their applications are also widespread. Detergents and mammalian respiration are two common examples in 
which surfactants are present. Living cell membranes are complex macromolecular assemblies comprised in 
large part of self-assembled phospholipids, of an amphiphilic nature 0. Sponge mesophases are formed as a 
result of an amphiphile dispersion or melt at an appropriate composition, and enjoy numerous applications in 
medical research as well as the pharmaceutical, cosmetic, food, and agro- and petrochemical industries 
Liquid-crystalline bicontinuous cubic mesophases of monoglycerides and the lipid extract from archaebacterium 
Sulfolobus solfataricus have been found at physiological conditions in cell organelles and physiological transient 
processes such as membrane budding, cell permeation and the digestion of fats j^. Amphiphilic cubic mesophases 
can also be synthesised for important applications in membrane protein crystallisation, controlled drug release 
and biosensors 0|- These phases are termed mesophases not only because their intrinsic internal length 
scales range between those characteristic of molecular and hydrodynamic (or macroscopic) realms, but also 
their mechanical properties are half-way between those found in a liquid and a solid P, 0, ■ 

Amphiphiles have the property of lowering the interfacial tension in a binary immiscible, say oil-water, 
fluid 8J . Given the bipolar nature of their molecular structure, amphiphile adsorption at the oil- water interface 
is a process which is energetically favoured relative to their entropically beneficial dispersion in the bulk. This 
effectively reduces the pressure tensor at the interface, making the immiscible species more alike. As more 
interfacial surface is created, so more amphiphile dispersed in the bulk can be accommodated at it. 

The effect of adding surfactant above a critical concentration to an oil-water mixture undergoing phase 
separation is to slow down the demixing process, which, with the addition of sufficient amphiphile, can be totally 
arrested. Langevin, molecular dynamics and lattice gas simulations have shown that, as the concentration of 
dispersed surfactant increases, the temporal growth law of the average size of the immiscible oil- water domains, 
of the power-law form t° in the surfactantless case is seen to cross over to a slower, logarithmic growth of 

the form (Ini)^, where a and 9 are fitting parameters and t is the time Qj^SE'l- Emerton et al. showed that 
increasing the surfactant concentration even further leads to growth well described by the stretched exponential 
A— B ex.p{—Ct^), where A, B, C and D are fitting parameters, including halted segregation at sufficiently late 
times 13]. Depending on temperature, pressure and fluid composition in such a stretched exponential regime, 
the amphiphile can self-assemble and force the oil-water mixture into a wealth of equilibrium structures. The 
self- assembling process is dictated by the competing attraction-repulsion mechanisms present among the species. 
Lamellae and hexagonally-packed cylinders are examples of these mesophases, also referred to as La and H, 
respectively, with continuous translational symmetry along one or two directions. Other examples are the 
sponge (L3) mesophase and the micellar (Q^^'^ or Pm3n, and Q^^"^ or Fd3m), primitive ("P", Q^^^ or ImZm), 
diamond ("D", "F", Q^^** or PnSrn) and gyroid ("G", Q^'^" or Ia3d) cubic mesophases, all of which lack 
continuous translational symmetry [lj|. Among all the aforementioned phases, only the sponge mesophase is 
devoid of long-range order and so cannot be classified as a liquid crystal: it is rather characterised by glassy 
features. 

A sponge mesophase formed by the amphiphilic stabilisation of a phase-segregating binary fluid mixture is 
called a microemulsion. Since we shall be dealing with oil and water in equal proportions, we shall be concerned 
with bicontinuous microemulsions. A bicontinuous microemulsion is a structure consisting of two percolating, 
interpenetrating oil and water phases separated by a monolayer of surfactant molecules adsorbed at the interface. 
Oil and water are isotropically mixed, and ordering is short range. Sponge phases formed by the dispersion of 
amphiphile in a single phase solvent differ from microemulsions in that it is a surfactant bilayer which underlies 
the structure, and the regions it divides are occupied by the same fluid component. A gyroid phase is also 
a bicontinuous, interpenetrating structure; however, ordering is evidently long range, whence its classification 
as a liquid crystal. In the gyroid, the locus where most of the surfactant resides is a triply periodic minimal 
surface (TPMS) whose unit cell is of cubic symmetry. The surface has zero mean curvature, no two points on it 
are connected by a straight segment, and no reflexion symmetries are present. Isosurfaces of the gyroid phase 
for which oil and water are not at equal composition (minority phases) form mutually percolating, three-fold 
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coordinated, regular lattices. Other examples of triply periodic surfaces of zero mean curvature arise in the 
P and D mesophases, the minority phase isosurfaces of which exhibit coordination numbers of six and four, 
respectively. 

The purpose of the present paper is to report on a theoretical study of the segregation kinetics in ternary 
amphiphilic fluids and the self-assembly of the sponge and gyroid mesophases. By progressively adding surfac- 
tant to an initially homogeneous immiscible oil-water mixture on the way to achieving arrested domain growth, 
we shall give an account of how the segregation kinetics of the fluid domains is affected by the addition of sur- 
factant, and study the features of the associated mesophases that are formed. The mesophases corresponding 
to such late time, arrested growth regimes are sponges which turn into gyroids as we increase the surfactant 
concentration. We shall also see that these phases exhibit temporal oscillations in the size of the oil-water 
domains, which we ascribe to Marangoni flows. 

II. OVERVIEW OF MODELLING AND SIMULATION OF AMPHIPHILIC FLUID SELF- 
ASSEMBLY 

Various methods have been used to date to model and simulate ternary amphiphilic mixtures and to study 
their phase segregation kinetics and the formation of microemulsions and liquid crystalline phases. We briefly 
review them in this section. 

Kawakatsu et al. studied segregation kinetics employing a two dimensional hybrid model with thermal noise 
but without hydrodynamics, combining a continuum, Langevin diffusion equation for the oil-water dynamics 
and Newtonian dynamics with dissipation for bipolar particles modelling the surfactant jlll |. They used a free 
energy in the form of a (/)^-Ginzburg-Landau expansion jisj plus terms modelling the surfactant-interface and 
surfactant-surfactant interactions. They found the average domain size of symmetric binary immiscible fluids 
with amphiphile to grow with time more slowly than t^/^, the latter expected for binary alloys in two and 
three dimensions. Laradji et al., instead of modelling the amphiphile as a particulate species, regarded it as a 
continuous density coupled to the oil- water order parameter in a 0'*-Ginzburg-Landau free energy p^ . In their 
work, they studied several cases of two-dimensional Langevin diffusion equations, one of which being the so-called 
Model D Model D incorporates noise, a conserved order parameter and surfactant density, but excludes 
hydrodynamics. Laradji et al. not only found logarithmic growth for the behaviour of the average domain size 
with time, but also observed a slowdown from it for higher surfactant concentrations and dynamical scaling 
for the structure function at intermediate times. Yao & Laradji, using a modified Lifshitz-Slyozov nucleation 
theory for continuum fields in two and three dimensions, studied how the Ostwald ripening dynamics of an 
asymmetric mixture of oil and water is altered by the presence of a surfactant species QJ^ . They found results 
similar to those of Laradji et al. [l^ . 

The segregation kinetics of amphiphilic fluids have also been studied with fully particulate methods such as 
classical molecular dynamics and, more recently, hydrodynamic lattice gases. Using a minimalist molecular 
dynamics model in two dimensions, Laradji et al. [ig found a crossover scaling function similar to previous 
Langevin and Lifshitz-Slyozov models [T^ . yet with a different algebraic exponent, and a slowing down 
from the algebraic growth laws for binary mixtures. Using two-dimensional hydrodynamic lattice gas models 
for symmetric 12i and asymmetric mixtures [l9| , the group of Coveney found that surfactant induces a crossover 
to a logarithmic slow growth, and, with sufficient surfactant, full arrest of domain growth which is well described 
by a stretched exponential function. The group found similar results with a three-dimensional hydrodynamic 
lattice gas model 20] . 



Particulate methods have also been used to tackle mesophase self-assembly. Using classical molecular dy- 
namics methods, Marrink et al. simulated evolution of a surfactant bilayer, initially set up on the morphology 
of a "D" TPMS, to study both the surfactant packing structure and how close such a bilayer would remain to 
the TPMS after relaxation Q. They, however, did not address self-assembly dynamics: time scales required 
for that are orders of magnitude above those reachable with atomistic techniques on present-day cutting-edge 
supercomputers. In dissipative particle dynamics (DPD) approaches, a Langevin dynamics with momentum 
conservation is solved to model ill-defined, mesoscopic dissipative particles interacting via repulsive, soft poten- 
tials; hydrodynamics is emergent and the amphiphile is represented by dissipative particles bound together by 
rods or springs 0i|2^|23. The DPD simulations of Groot & Madden of copolymer melts showed that melts 
of symmetric amphiphile led to lamellar phases, whereas a gyroid-like structure appeared only for asymmetric 
amphiphile as a transient phase, precursor to a hexagonally packed tubular phase. Nekovee & Coveney, using 
the lattice-Boltzmann model we employ in this work, were able to reproduce the "P" mesophase in a binary 
amphiphilic mixture of surfactant and solvent |24j . 
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Many of the simulation studies on the formation kinetics of microemulsion and liquid crystalline mesophases 
have made use of stochastic Langevin diffusion methods, in which mass currents are driven by chemical potential 
gradients computed from free energies of the ubiquitous ^''-Ginzburg-Landau expansion form. These models 
treat the amphiphile only implicitly through the functional dependence of the surface tension parameter with the 
amphiphile density 25, 26, 27, 28, 29J. In cases in which the amphiphile is a copolymer, however, the free energy 
is often derived from polymer models which aim at accounting for the amphiphile's molecular structure with 
a certain degree of specificity §3, 1^ ■ The validity of these Flory-Huggins type approaches rests on being 
able to derive the free energy from a microscopic model of the complex fluid mixture, which not only might 
entail considerable difficulty but does require the segregation to be a quasi-static, local equilibrium process. 
Under general far from equilibrium conditions, such as occurs in the sudden-quench scenario so often employed 
in the literature, equilibrium thermodynamic potentials are known not to adequately describe the process. 
Besides, free energy based methods also require surfactant adsorption and relaxation on the interface to be 
much faster than interface motion, a so-called adiabatic approximation. Free energy approaches are frequently 
represented as paradigms of thermodynamically consistent mesoscopic methods; some of them also pursue 
chemical specificity in elaborate empirical exercises amounting to little more than parameter fitting of polymer 
models. The philosophy behind them, nonetheless, is the use of macroscopic, local equilibrium information to 
specify a stochastic, and hence mesoscopic, non-equilibrium dynamics. None of these methods offers a dynamics 
satisfying detailed balance, let alone an i/-theorem (Lyapunov function) guaranteeing irreversible evolution 
towards the equilibrium state described by the prescribed free energy. As a consequence, the 'thermodynamic 
consistency' of these methods remains on shaky CTOunds. 

The fact that some free energy approaches [3ll l32i | focus on the specific molecular structure of the amphiphile 
raises the question of what use particulate methods, such as is the one we report in this paper, have in the 
simulation of amphiphilic fluid systems. Our method, by reducing the description of the amphiphilic molecule 
to its minimal possible expression — a point dipole, retains the minimum number of degrees of freedom necessary 
to model interfacial adsorption and micellisation, and, additionally, in a hydrodynamically consistent framework 
which does not require processes to be quasi-static. With these basic properties at our disposal, we want to fully 
exploit our model's capabilities to determine the non-equilibrium amphiphilic dynamics and the equilibrium fluid 
structures arising from it. The minimalistic bottom-up approach is in line with the fact that, far enough from 
criticality, distinct molecular structures and microscopic dynamics can produce similar macroscopic behaviour — 
this is universality emerging from microscopic complexity ^33,]. In addition, particulate methods are much more 
suitable than conventional continuum fluid dynamics methods |34j for the simulation of interface dynamics. 
Such dynamics is an emergent property of the underlying inter-particle interactions among the immiscible 
species; a set of continuum partial differential equations describing the locus of the interface is, rather, its 
macroscopic manifestation, and its solution a much more labourious endeavour. A fortiori, modelling surfactant 
adsorption and self-assembly in an explicit fashion via particulate methods provides a more realistic picture of 
the microscopies than doing so at the continuum, macroscopic limit described by free-energy approaches. 

Lattice-Boltzmann (LB) methods were originally developed as a means of reducing the computational cost 
associated with lattice-gas automaton (LGA) algorithms 35]. LB methods evolve a single-particle distribution 
function via a discretised Boltzmann equation, usually in the linearised, relaxation-time (BGK) approximation. 
Such a single-particle distribution, at a particular time slice and spatial position on the lattice, is an average over 
the LGA velocity space for a statistically large number of different microscopic realisations (initial conditions). 
The fact that much of the phenomenology of binary immiscible and ternary amphiphilic fluids occurs for 
small spatio-temporal gradients permits us to take the mean-field (or molecular chaos) approximation, and the 
Boltzmann-Grad limit in which such an approximation holds, as heuristically appropriate for the modelling of 
their universal properties. Heuristics come into play in that tunable parameters are introduced in LB models 
in order to reproduce desired quantities of dense and/or complex fluids, such as surface tension, viscosity 
and thermal conductivity (for required values of Reynolds and Prandtl numbers), stress tensors (for required 
viscous or viscoelastic behaviour), and equations of state (for liquid-gas and phase segregating transitions). It 
is worth noting that the increasing popularity of LB methods in recent years is primarily based on pragmatic 
considerations associated with their simplicity and algorithmic efficiency. 

This paper presents the first quantitative account of amphiphilic phase segregation dynamics using a three- 
dimensional model based on the Boltzmann transport equation. It describes the spontaneous self-assembly 
of the gyroid liquid crystalline cubic mesophase and an order-disorder transition between the latter and the 
sponge mesophase, of glassy features. The remainder of the paper is structured as follows. In Section III we 
describe the model. In Section IV we look at how the segregation kinetics of the fluid domains is affected 
by the addition of surfactant. Section V studies the temporal oscillations of the average domain size and 
the structure function, which are only observed for segregation-halted regimes. In Section VI we characterise 
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the morphology of the mesophases corresponding to those regimes via direct- and Fourier-space imaging, and 
identify the sponge<->gyroid structural transition. Finally, we provide conclusions in Section VII. 



III. A LATTICE-BOLTZMANN MODEL FOR TERNARY AMPHIPHILIC FLUIDS 



The amphiphilic lattice-Boltzmann model we employ in this paper is derived from that originally proposed 
by Chen et al. 36, 37]. The method can be regarded as a fully mesoscopic, bottom-up approach, which does 
not require the existence of a thermodynamic potential describing phase transitions. In fact, the method 
is athermal in the sense that, for algorithmic efficiency reasons, the microdynamics is devised ex professo 
to conserve velocity moments of the distribution function only up to first-order; this simplification is valid 
wherever thermal fluctuations are negligible, e.g. away from criticality. This is, for example, the case of deep 
quenches into the spinodal region of the fluid's phase diagram, which is our case in this paper. As opposed 
to top-down LB methods, based on the imposition of a free energy functional H^H^, the global dynamics 
arise as an emergent property of the interactions between mesoscopic levels of description, in agreement with 
a complexity paradigm |33| . Oil- water segregation is achieved via inter-species forces which modify the fluid's 
macroscopic velocity. The dynamics in the bulk of each binary immiscible species (e.g. oil and water) can be 
derived from a Boltzmann equation with a forcing term. An amphiphilic molecule is modelled as a continuously 
orientable massive point dipole subjected to thermal noise and relaxing towards an equilibrium that minimises 
its interaction energy with mean fields generated by its nearest neighbours on the lattice. The densities of 
surfactant, oil and water evolve via coupled lattice-BGK equations. This is a mean-field approach which exhibits 
Galilean invariance and reproduces correct hydrodynamics. We have also shown, in a previous paper which 
serves as a reference benchmark for this study 0, that the model reproduces the dynamical scaling hypothesis 
during the phase segregation experienced by binary immiscible (oil-water) fluids. Its algorithmic simplicity 
allows it to achieve extremely high performance on massively parallel computers j4Cl( | , and substantially reduces 
the domain of numerical instability present in free-energy-based LB methods |9|. Because an if-theorem is 
lacking in essentially all multiphase lattice-Boltzmann models hitherto proposed we consider it artificial 
to try to enforce a prescribed thermodynamic equilibrium in these schemes; a method which is algorithmically 
simpler, fully mesoscopic, mean-field and bottom-up is of greater fundamental interest. 



A. Binary immiscible fluids 



The core of our model is a lattice-BGK equation governing the evolution of the mass density distribution 
TO"n^(x, t) of component a in an interacting fluid mixture at position x, instant t, and for discrete molecular 
velocity Cfe, on a regular lattice and in discrete time. Here, m" is the particle mass which we set to unity for con- 
venience, and the single-particle distribution n^(x, i) obeys the lattice-BGK relaxation-streaming mechanism: 



<(x + Cfe,t+l)-<(x,t) =17^ (1) 

where the collision term has two contributions accounting for the kinetics of non-interacting (ideal) plus inter- 
acting (non-ideal) multicomponent species, respectively: 

17?(x, t) ^ f7i")"(x, t) + E E < ■ (2) 

a I 

the sums extending over all available species and directions, and 

17W"(x,t)^-"^^^'^^~5^"''("'^^ (3) 

Here, the time increment and lattice spacing are both unity, x is a node of such a lattice, a = r, b (e.g. oil (r) 
or water (b)), and is one of the 24 (= A'^vec) discrete velocity vectors plus one null velocity of the projected 
face centred hypercubic D4Q25 lattice we use to guarantee isotropy in the macroscopic equations that the 
model reproduces for a bulk, single phase fluid The parameter r" defines a single relaxation rate towards 
equilibrium for component a, A^" can be regarded as a matrix element of a cross-collision operator A which is a 
function of both r" and the acceleration a", the latter being experienced by a fluid element due to its neighbours. 
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as will be defined later. The function n'j^^'^'^\yi, t) in Eq. (0 is the discretisation of a third-order expansion in 
Mach number of a local Maxwellian representing the local equilibrium state of the ath component, 



Q(oq) 



(x,t) =tjfen°(x,i) 



1 



rCfc • u 



2c! 



2c2 



6c? 



(Cfc • u) 



2c4 



U^(c/c • U) 



(4) 



where LUk are the coefficients resulting from the velocity space discretisation, and Cs is the speed of sound, 
both of which are determined by the choice of the lattice. For the projected D4Q25 lattice we use, the speed 
of sound is Cs = l/VS, oJk = 1/3 for the speed Ck = 0, and 1/36 for speeds Cfc = 1 and -^72. In Eq. 
u = u(x, t) is the macroscopic velocity of the mixture, through which the collision term couples the different 
molecular velocities c^. This is because u is a function of the components' macroscopic velocities, defined as 
n"(x,t)u" = ^,<(x,t)cfe. 

A judicious choice of the coefficients in the expansion of the equilibrium distribution Eq. allows for mass 
and momentum to be (non-locally) conserved for the non-interacting, ideal gas mixture case, i.e. 



(0)q _ 



0. 



(5) 



It can be shown that in the limit of creeping flows to second order, i.e. « 0, the expression for the fluid 
mixture's macroscopic velocity u required for momentum conservation in the absence of interactions, as a 
function of u" , simplifies to that obtained for a second-order expansion of the equilibrium distribution, namely 
u = /X^a f^' which we have incorporated in our implementation. 

The form of the collision term ^ derives from adding an increment Au" to the fluid mixture's macroscopic 



velocity u which enters in the equilibrium distribution i.e. ^^^(u) 
and a" = F"/p". Here 



r2^°^"(u + Au") where Au° = a"r° 



F"^^(x, t) EE -ri^, t) E ffaa E ^''(^'' *)(^' - ^) 



(6) 



is the mean- field force density felt by phase a at site x and time t due to its surroundings; gaa is a coupling 
matrix controlling the interfacial tension between the fluid species, interface adsorption/desorption properties 
of the surfactant molecules, and the surfactant-surfactant interaction; tp"' is an effective mass which serves as 
a functional parameter and can have a variety of forms for modelling various types of fluids. We only allow 



n"(x, t) , where = ^ 



nearest- neighbour interactions, x' = x -I- Cfc, and choose ■ip°'{'x.,t) = 1 — exp 

This choice for has also been made by Shan and Chen to model liquid-gas phase transitions [4^ 
we shall see, our motivation here is different. 



fc "fc- 



although, as 



B. Amphiphilic fluids 

The incorporation of a third, amphiphilic species not only requires the inclusion of an extra variable ("s") 
for the superscript denoting the species in Eq. but also a modification of the cross-collision operator A 
since amphiphiles interact with fluid elements and between themselves. In addition, the physics of amphiphilic 
molecules, namely, self-assembly and adsorption to immiscible fluid interfaces, cannot be modelled without 
introducing a new type of body force: in Subsection IIII Al ordinarv bulk fluid species are thought of as point- 
like particles given that their interactions depend on their relative distances alone. For surfactant molecules, 
however, their orientations are important too |36| . and a dipole is the simplest configuration to mimic their 
essential character. In short, we must extend the scalar lattice-BGK model hitherto described into a vector 
model. 

Each surfactant molecule is represented by an average dipole vector d(x, t) at each site and time step, whose 
orientation is allowed to vary continuously. The average is taken over nearest neighbours before advection, 
according to the propagation equation 

n^(x,i + l)d(x,i+l) = E^fc(x-Cfe,t)d(x-Cfe,i), (7) 

fc 

where the tildes denote post-collisional values, as defined by Eq. ^ for the A = {gas = 0) case if we replace 
the leftmost summand with h'^ (x, t) . For the sake of simplicity and computational efRciency, the model does 
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not assign microscopic velocities to single dipole vectors but to site-averaged surfactant densities instead, as 
can be seen, for example on the right hand side of Eq. 0. 
Dipole relaxation is governed by the BGK process 



d(x,i) = d(x,i) 



d(x, t) - d'=i(x, t) 



(8) 



where is a new parameter controlling the relaxation towards the local equilibrium d°'^(x, t), which is under- 
stood as the average orientation with respect to the Gibbs measure, i.e. 



(9) 



where d^J7 is an element of solid angle whose director is the unit vector Q. representing the dipole orientation, 
and P is an inverse temperature-like parameter. The modulus of the distribution ^ ranges between and the 
scale value do, chosen to be unity for convenience. That, along with > 1, guarantee the magnitude of the 
dipole vector to be less than do at all times. Equation favours surfactant orientations which minimise the 
energy = —tt ■ h(x, t), where h(x, i) is the sum of the mean fields created by surrounding bulk fluid and 
surfactant, namely 



h'=(x, i) = ^fc^n"(x + Cfe,t)cfc , 

a k 

h^(x,i) = ^ [n^fe(x,i)d(x,i) +^7i^(x + cz,i)e; •d(x + Q,<) 



(10) 
(11) 



allowing for nearest-neighbour interactions only. The first equation is a discrete approximation to the colour 
gradient for the immiscible species, where qa = 0, ±1 is the colour charge of species a. The second equation 
is a dipole vector density, where summation over k performs local dipole averaging, summation over I includes 
all nearest neighbour contributions, and the second-rank tensor 0; = I — DciCijc^ , where c is the modulus of c; 
and D is the spatial dimension, picks up desired orientations from nearest-neighbour dipoles. Finally, Eq. 



do 



coth(/3/i) - 



1 



h, where h is the magnitude 



can be integrated analytically in three dimensions to give d'^'^ 
of h and h its unit vector. 

The new interactions that modify the interspecies collision operator A are the force on an immiscible fluid 



element from other fluid elements and amphiphiles, F" = 5brF° 



where F"''^ is that in Eq. (jSJl, and 



the force on an amphiphilic molecule from neighbouring fluid elements and amphiphiles, F'^ = ^bsF®''^ + gssP^'^- 
In these expressions, ffbs and g^s are coupling scalar parameters, and the analytical expressions for each 
force term, derived in [3a, are 



-2.9„,^"(x,t)^d( X + ci,t)Oi ■ V/(x + ci,t) 



2i;%K, i)d(x-HCi,<) 

4D , „ , . 
2-5ssV' (x, i) 2^ 



9as fi'/c -0" (x Cfc , t) 



|d(x + Cfc,i) - Ok •d(x, t)ck 



d(x -t- Cfc, t)d(x, t) + d(x, i)d(x -I- c^, t)] •c,}^^(x + Ck,t) 



(12) 
(13) 

(14) 



Equations (|12|l . (|13|l . and (|14|) were derived considering only nearest- neighbour interactions, modelling each 
dipole as a dumbell of oppositely colour-charged particles displaced ±d/2 from the dipole's centre of mass 
location x, and carrying out Taylor expansions of the force lO to leading order in d about x as well as those 
at the neighbouring sites [s^. Also, Eq. H13|) is the reaction to force H12|l . and Taylor expansions in the 
ratio of |cfc| to the length scale that the colour gradient sets can be used to further simplify the expressions. 
Finally, additional coupling parameters gaa have been introduced, where gss should be chosen negative to model 
attraction between two amphiphile heads or tails, and repulsion between a head and a tail. 
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C. Selection of the parameters for the simulations 



The model is implemented as a parallel code in Fortran90 making use of the Message Passing Interface 
parallel paradigm |44| and spatial domain decomposition, and incorporating wrap-around, periodical boundary 
conditions. It was executed on 16 to 64 processors on SGI Origin2000 and Origin3800 parallel platforms. The 
form -0 = 1 — exp[— n(x, t)] for the effective mass in the force in Eq. © was heuristically chosen so as to broaden 
the region of numerical stability in parameter space: numerical instabilities can arise as the result of high values 
of forces and speeds, and are more likely to occur in our model when surfactant interactions are included than 
for binary immiscible fluids 0, • 

Preliminary studies allowed us to determine the values of the model's various parameters for which an initially 
thorough mixture of two immiscible fluid phases plus a dispersed amphiphilic species produced a segregated 
mixture with arrested domain growth |45| |. Those values were the surfactant thermal parameter (3 — 10.0, 
all particle masses and relaxation times set to 1.0, and coupling constants ghr = 0.08, g^s = —0.006, and 
gss = —0.003. (Masses, m", enter in the description through p°'{x,t) = ^k'm"'n'^{y^,t)-) We simulated the 
behaviour of a ternary mixture by varying the coupling constants around the values mentioned above, and for 
initial surfactant particle densities ranging in the interval 0.00 < rS^^^ < 0.90. The lattice sites and directions 
were initially populated with flatly distributed mass densities, < p^(x, 0) < m"n(°)"/A^vec, where nf^"^" is 
the particle density of phase a, and k numbers each of the A'vcc = 24 velocity vectors. In addition, we used 
periodic boundary conditions in all three dimensions. We determined that setting n^^^" > 0.6 for both species 
guaranteed immiscibility. In all the simulations we present here we set oil: water mass fractions to 1:1, specifically 
at n(°)" = 0.7 for a = r,b. 

It is experimentally known that the addition of amphiphile into an immiscible fluid mixture reduces the 
interfacial tension, as has also been reported for various lattice gas models in two and three dimensions [l3Ll20| . 
To confirm that our model reproduces this important property, we ran simulations on a 4 x 4 x 128 lattice of a 
planar interface with surfactant adsorbed onto it and whose initial density was varied between simulation runs. 
The surface tension was calculated with the line integral along the normal to the interface 



cr = 



+ 00 

' PM ~ P.,,{z) dz , (15) 



where, for the pressure tensor P = {Pij}, we used the expression 

a k 
Q,a x' 

We restrict ourselves in this study to nearest-neighbour interactions, x' = x + c^, and transversal symmetry 
allows the second summand within the integrand in Eq. (|15|l . which in general is ^^iPxx+Pyy), to be simplified as 
shown. Equation (|16|l contains a kinetic (first) term, the momentum flux, due to the free streaming of particles 
corresponding to an ideal gas contribution, plus a potential or virial (second) term due to the inter-particle 
momentum transfer derived from the force ® 47| E^l . 

Figure^ shows the surface tension a plotted against initial surfactant density, and details on parameters and 
densities used are included in the caption. Notice that in the regime the binary fluid is in, and for the values 
of surfactant density we use, the surface tension decreases linearly with surfactant concentration. It is entirely 
possible that there may be departures from linearity were we to increase the surfactant concentration beyond 
that shown in Fig. ^ because of interfacial saturation with surfactant, as observed in two and three-dimensional 
lattice-gas studies ^ . 



IV. DOMAIN GROWTH KINETICS 



We ran simulations starting with a homogeneous mixture of oil and water particles mixture to which surfactant 
was randomly added across on the lattice. Lattice sizes employed were 64^ and 128'^ to assess finite size 
effects. Each lattice site was populated with a density uniformly distributed in the range zero up to the values 
summarised in Table 
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The average size of the oil-water domains is a natural measure of the degree of segregation within the 
mixture. We define it as the inverse first moment of the spherically averaged oil-water structure function, 
L{t) = 2'K/ki{t), where fci(t) = kS{k,t)/^f, S{k,t). The spherically averaged oil-water structure function, 
S{k, t), is S'(k, t) I 1, where the ^(k, t) is the oil-water structure function, 

S{)^,t)^\-<^^{t)\ (17) 



V 



and denotes summation over the set of wavevectors contained in the spherical shell n — < k < n -I- ^, 
for integer n. Equation (|17|l is the Fourier transform of the spatial auto-correlation function for the oil-water 
order parameter = — p° , where V is the lattice volume, the volume of the lattice unit cell, and (l)'^{t) 
is the Fourier transform of the fluctuations of the order parameter, 0. Our choice of the structure function, 
rather than alternative measures of domain size such as the auto-correlation function, was made on the basis 
that it is directly proportional to X-ray or neutron scattering intensities, hence facilitating direct comparison 
with empirical data [49j. 

In Fig. El we plot the temporal evolution of the average domain size L for the surfactant concentrations of 
Table The amount of surfactant needed to slow down the kinetics of the binary immiscible oil- water mixture 
(simulation run 01) is seen to be relatively low. We now need to find the growth laws that best fit these data. 
Previous simulation studies, for 1:1 oil- water fluid mixtures with or without surfactant 0,0,113,133, have 
found algebraic, logarithmically slow and stretched exponential behaviours, as follows 

ai{t-hiY\ (18) 
a2(lnt)^= , (19) 

ag-fesexp [-C3(t-d3)"='] , (20) 

to be those characterising the temporal growth of the average domain size, L{t)^ of an oil- water mixture without 
surfactant, Eq. 118|l 9], and when surfactant is added above a minimum threshold concentration, Eq. H19|l . and 
at a sufficiently high amphiphile concentration, Eq. H2U|) , the latter being a regime for which arrested growth is 
reached at late times. The coefficients ai and hi [i — 1,2,3) are fitting parameters. While we shall take these 
functional forms as suggested choices, we would also like to find out how closely they in fact fit our data. 

Linearity in the {t, L) data cloud on a log-log plot would permit us to ascertain whether or not the data 
follow Eq. H18|l . regardless of the zero-time offset value bi since this is a horizontal displacement. To find out 
which data may be better fit by Eq. (O, we would require (logi, L) pairs of data in a search for linearity on a 
log-log plot. This method, however, is not likely to be of much help given the small difference between plots of 
the logarithm of a data series and plots of the logarithm of such logarithmic data, as we shall see. We therefore 
prefer to adopt the criterion of considering candidates for the model of Eq. H18|) from the log-log linearity 
method, while resorting to both visual inspection and a search for a reduced chi-squared statistic (x^/ndf) close 
to 1.0 in order to identify a slower growth such as that of Eq. (|19|l (ndf is the number of degrees of freedom). 
Finally, Eq. H2U|) possesses a distinctive horizontal asymptote which best fits data whose domain growth at late 
times is fully arrested. 

From the linearity of curves in Fig. |31we can infer that simulation runs 01 and 02 (Table|l|) are good candidates 
for the growth model of Eq. (|18l) . Figure O however, leads to the same conclusion, as expected given the small 
difference between these two plots. We then resort to looking at the x^/ndf statistic in assessing how well Eqs. 
(|18|l and H19|l fit simulation runs 01, 02 and 03, see Table HH The binary immiscible fiuid simulation run 01, 
with no surfactant present, exhibits an exponent consistent with the system being in a crossover between the 



simulation run 


01 


02 


03 


04 


05 


06 


07 


08 




0.0 


0.15 


0.22 


0.30 


0.35 


0.40 


0.60 


0.90 


a;" 


0.0 


0.21 


0.31 


0.43 


0.50 


0.57 


0.86 


1.3 



TABLE I: Surfactant densities employed in the study of the algebraic-to-logarithmic and logarithmic-to-stretched expo- 
nential transitions. The mass fraction, x^, is the ratio of n^-^^'' to n'"^'' = n^"'"^ = 0.7, and the rest of the parameters used 
were g^r ~ 0.08, Qhs ~ —0.006, gss ~ —0.003, masses and relaxation times set to 1.0, and /3 — 10.0. The lattice used was 
sized 64"^ for all simulation runs except 01, for which it was 128'^ in order to avoid finite size effects entering at about 
L « 25. 
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known diffusive {t^^^) and viscous hydrodynamic (t^") regimes, already reported for binary immiscible fluids 
simulated with the lattice-BGK model we employ in this paper 9] . Simulation run 02 has the peculiarity that 
Eq. p8|l holds (poorly) only during an initial transient, and Eq. H19|l takes over to give a very good fit at later 
times, t > 1100. This transient is due to the time required by the surfactant to adsorb onto the interface and 
affect the binary immiscible interfacial dynamics. Finally, simulation run 03 is best fit by Eq. (UD), although 
the high x^/ndf value indicates that the data contain more detail than the model does. In addition, from 
Fig. this mixture segregates at a slower speed than that given by Eq. H19(l . yet it does not reach total 
arrest, at least up to 7200 time steps. Rather, total arrest is seen at higher surfactant concentrations, as in runs 
06 and 07 (see Fig. 0). We conclude that simulation run 03 represents a fluid which is in a transition regime 
between the logarithmic and the stretched exponential behaviours. A similar behaviour was previously observed 
by others using lattice-gas methods in two ^3 and three-dimensions ^^'^ lattice-Boltzmann methods in 
two dimensions 37] . Emerton et al., using a two-dimensional lattice-gas model, reported the divergence of the 
coefficients of Eq. (|20|l in an attempt to fit data for which total growth arrest had not been achieved ^3 ■ The 
fits to our data, which include error bars, also produced the same divergences. Their fluid mixtures as well as 
ours, we conclude, were, rather, in a transitional regime well described by a growth law slower than Eq. (|19|l 
which still allowed for domain growth. It is, however, possible that growth arrest could be achieved at later 
times; this pre-arrest regime would then be a long-lived transient. 



simulation run 


Cl 


xVndf 


C2 


xVndf 


01 


0.896 ± 0.007 


0.18 






02 


0.644 ± 0.004 


7.5 


3.850 ±0.010 


0.92 


03 






2.649 ±0.022 


39 



TABLE II; Fits of the average domain size growth with time to the models of Eqs. 11811 and 11911 for simulation runs 
corresponding to surfactant mass fractions 0.0, 0.21 and 0.31, from top to bottom, respectively, as detailed in Table U 
Lattice sizes used were 128^ for simulation run 01 and 64"^ for the rest. Poor fits are indicated as blank fields. Simulation 
run 02 shows two behaviours in its temporal evolution, Eq. 1181 1 for t < 1100 and Eq. 11911 for t > 1100. Note the very 
good value of the x^/ndf statistic for the latter. The poor value of the statistic for simulation run 03 indicates that 
Eq. 11911 is insufficient and a more detailed model is required, albeit not Eq. 12011 . 

We now look at wavenumbers of the spherically averaged structure function, S'(fc), other than the first moment, 
already provided by L{t). Figure 0] shows the spherically averaged structure function for simulation run 06 at 
several time steps. The temporal evolution of the curves 

resembles the segregation kinetics for binary immiscible fluid mixtures, except that domain growth arrest for 
late times makes them tend to superimpose. Note that a hump appears at these times, indicating the formation 
of structures, statistically weak, of size close to half the lattice side length. Inspection of 0(x) snapshots suggests 
the spurious presence of elongated domains of such sizes which are extended rather than folded. At the late 
times we examined, these elongations tend to vanish or fold. Still in Fig. ^ it is worth noting that for all 
length scales above a threshold (about k < 0.9), curve superposition is not sharp. This is a consequence of 
the fact that for the fluid composition of simulation run 06, and those of higher surfactant concentrations, 
there are small temporal oscillations in S{k). All these mixtures have in common that they have achieved total 
growth arrest — in fact, L{t) decreases in time for simulation run 08, as we shall see later on and discuss in more 
detail. Oscillations in the structure function and a hump at low wavenumbers have been reported previously in 
a hydrodynamic Langevin model of sponge phase dynamics, using field-theoretic methods |2S| |. However, this 
approach did not consider the amphiphile concentration explicitly but, rather, embedded it into a Ginzburg- 
Landau free energy through the surface tension, in a scenario where amphiphile relaxation is assumed to be fast 
compared to that of the oil-water order parameter. 

In Figs. andlSb we show the spherically averaged structure function at time step 7200 of the mixtures in 
Table m for the larger and smaller length scales, respectively. As the initial density of surfactant is increased 
in a series of replica initially homogeneous water-oil-surfactant fluid mixtures, as indicated by Table ^ it is 
expected that the oil-water structure function peaks will move to higher wavenumbers, decrease in intensity, 
and broaden 0,11, H^l- This is indeed what we observe in Fig. Note that at smaller length scales. Fig. [SJa, 
the exponential decay of the structure function that occurs for simulation run 01 does not hold for the ternary 
amphiphilic mixtures. This can be explained by the contribution of small micellar structures that form in the 
bulk of each immiscible phase, more likely to take place for mixtures of higher surfactant concentration. Indeed, 
in Fig.lSja, the latter exhibit the most manifest deviations. 
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V. SELF-SUSTAINED OSCILLATIONS 



Arrest of domain growth occurs for high surfactant density only, cf. Fig. [51 as expected. Further inspection, 
however, shows that not only are there small temporal oscillations of the average domain size, as we mentioned 
at the end of the last section, but also that they do not die out during the simulation window. Similarly 
to what was previously reported using a bottom-up lattice-Boltzmann method in two dimensions akin to the 
one employed here j37|. the amplitude of the oscillations is very small compared to the average domain size, 
and smaller than previous lattice-gas simulations in two and three dimensions l2(il |. The fact that, in 
lattice-gas methods, these oscillations persist after ensemble averaging is consistent with their occurrence in 
lattice-Boltzmann approaches, since the latter are effectively ensemble-averaged versions of the former. Since 
the systems we simulate are dissipative and isolated (there is no mass or momentum exchange with external 
sources), oscillations, however, are expected to die out at sufficiently late times. 

Motivated by the observation of oscillating average domain sizes, we performed additional simulation runs in 
order to check the role of the coupling constants gss and gbs in the reproduction of such oscillations, our hypoth- 
esis being that both an increased surfactant-surfactant interaction and an increased tendency for surfactant to 
adsorb on the interface might be expected to have an influence on their frequency and amplitude. In Table 
mil we summarise the parameters used in the new simulation runs (09 and 10) along with those of previous 
'oscillating fluid mixtures.' 



simulation run 




1= 


3ss 


3bs 


06 


0.40 


0.57 


-0.0030 


-0.006 


07 


0.60 


0.86 


-0.0030 


-0.006 


08 


0.90 


1.3 


-0.0030 


-0.006 


09 


0.60 


0.86 


-0.0045 


-0.006 


10 


0.60 


0.86 


-0.0030 


-0.009 



TABLE III: Parameters employed in studying domain size oscillations, whose onset occurs for surfactant mass fractions 
x" > 0.57; the remaining parameters of the model are stated in the caption of Table |I] also for the additional runs 09 
and 10. In lattice units. 

Figure |H| shows the temporal oscillations in the average domain size for the mixtures of Table IIIII and 
Fig. [7| shows their structure functions at time step 17000. All these mixtures exhibit domain growth arrest; 
interestingly, Fig. 1^1 shows that the average domain size shrinks in time for some of them (mixtures 08 and 
09). In addition, we uncover the role that the coupling constants gss and ghs have in the oscillations: whilst 
increasing j^ssl seems to enhance their frequency, an increase in \ghs\ drastically dampens them and reduces 
their amplitude. However damped the oscillations of simulation run 10 may seem, zooming into smaller scales 
reveals the existence of minute oscillations (less than O.IO lattice sites in amplitude), which is not the case 
for simulation runs 01 through to 05. (Note that the length scales reported in Fig. |Hlare lattice averages; an 
amplitude being less than one lattice site hence remains physically meaningful.) Oscillations are, therefore, the 
signature of all growth-halted regimes. 

The structure function plots of Fig. |7| provide further insight into the role of coupling constants gss and g^s 
in the oscillation dynamics. Note that mixture 08 produces a peak of intensity similar to that of mixture 07, a 
feature already seen at much earlier times (see Fig. |5^). This peak height similarity could have been ascribed 
to a transient, such as turned out to be the case for the difference in peak intensities between mixtures 06 
and 07; however, it persisted in time. Mixture 09 also shows a peak intensity similar to that of mixture 07. 
Peak intensities bear a direct relation to the steepness of oil-water domain walls and, hence, to their surface 
tension. The fact that increasing the surfactant concentration (in mixture 08 compared to mixture 07) does not 
reduce the surface tension denotes that the interface is close to its saturation limit with respect to surfactant 
adsorption. If enough surfactant is dispersed in the bulk, a process of diffusion towards and adsorption onto 
the interface could continue to occur, much slower compared to the initial adsorption leading to growth arrest, 
which could explain the slow domain size reduction. In the cases of simulation runs 08 and 09, close to interface 
saturation, surfactant concentration in the bulk is high. An amphiphilic mixture being close to the saturation 
limit implies that the value of its surface tension is the lowest among all amphiphilic mixtures sharing the same 
composition, relaxation times and coupling constants gas- Surface tension may be further reduced only by 
allowing more surfactant molecules onto the interface, which can be done by increasing |gbs|- This is exactly 
what we observe in Fig. for fluid composition 10. 
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As we saw in Fig. 0] small oscillations in the average domain size indicate that the structure function varies 
in time back and forth between distributions of sizes which are close to each other. The first moment of such 
distributions, as studied in Fig. El may not be representative of the dynamics at other length scales, as we shall 
see immediately. In Fig. |Slwe show the temporal evolution of S{k) for mixture 09 for a range of wavelengths. 
Note three characteristic features of the S{k) curves: they all oscillate, decrease for k < 0.785 and k > 1.08, and 
increase or remain stationary in the long time average for k « 0.884 and k « 0.982. This behaviour corresponds 
to the sharpening of the distribution S{k) with time. Modes with k > 1.28 {L < 4.91) decay fast enough 
< 0.1 for t « 1000) for them to be negligible in terms of their contribution to the fluid mesostructure. 
Other decreasing modes take much longer {t > 30000) to vanish. 

Our study of the oscillations would be incomplete without looking at frequency power spectra. The time 
series we analyse correspond to S{k — 0.589) of fluids 06, 07, 09 and 10; this choice is made on the basis that 
this wavenumber apprehends characteristic features of each data set. From each time series we subtracted its 
longest waves (i.e. its envelope), computed as the average j J2t' ^i^j^'), A being a lag large enough so as to 
decouple high-frequency from low-frequency waves (A — 5000 time steps), the sum extending over the interval 
t— A/2 < t' < i — A/2. The Fourier transform of the resulting time series we take as the definition of S(k, uj), cf. 
Fig. |51 Note therein two high peaks for simulation run 06, and a collection of weak peaks (which we define as 
those whose heights are less than 5% the height of the largest peak) occuring for higher frequencies. An increase 
in surfactant density (simulation run 07) causes the number of excited high-frequency modes to grow slightly, 
yet they also decrease in intensity. Simulation run 09, which differs from mixture 07 in having an increased 
\gss\i very clearly exhibits a substantial increment in the number of excited high-frequency modes. Finally, the 
spectrum for mixture 10 corroborates the quenching effect on fluctuations caused by increasing \ghs\- 

The term Marangoni instability describes a convective flow caused wherever an inhomogeneous temperature 
or mass distribution locally alters the interfacial tension |5lj | . By visualising the oil- water interface for mixtures 

06 to 10 we observed that the density of adsorbed surfactant is not evenly distributed on it; hence the conditions 
are set for the appearance of Marangoni instabilities. Figure EH displays the late time evolution of a subdomain 
of a fluid of the same composition as simulation run 09 but simulated on a larger (128'^) lattice. We display 
the surfactant density on a slice through the mid-plane of the subdomain, along with the locus of the oil-water 
interface depicted as an isosurface cropped close to the plane. Surfactant inhomogeneities on the interface are 
evident from these images, as well as the existence of a slow, creeping flow. Distinctive features include the 
regularity of the order parameter (which we shall study in detail shortly) , the existence of high surfactant density 
necks bridging adjacent portions of the interface, and local regions where regularity is absent, reminiscent of 
the defects in crystalline materials, which possess their own larger-scale dynamics. 

VI. THE SPONGE AND GYROID MESOPHASES 

It is known that, in an initially homogeneous 1:1 oil-water mixture, the arrest of phase segregation experi- 
enced through the addition of sufficient amphiphile can lead to the formation of a thermodynamically stable 
bicontinuous sponge phase 0,0,|23- In Fig. ^2 we show the late time morphologies for fluid compositions 06, 

07 and 08. They are displayed as the 0(x) = 0.37 isosurface, corresponding to a water-in-oil, "rod-like" scenario, 
where water is a minority phase and oil is in excess (the order parameter ranges as —0.69 < (j) < 0.68 over 
the lattice at that time slice) . The structure suggested by minority-phase isosurfaces and the structure of the 
interface ((/'(x) = 0) for fluid composition 06 resembles that of a microemulsion, for which structural disorder is 
the predominant feature. Fluid composition 08, by contrast, shows an evident resemblance to minority-phase 
images seen in transmission electron microtomography of the gyroid "G" cubic phase |52l| . The morphology is 
an interweaving, chirally symmetric, three-fold coordinated, bicontinuous lattice. Fluid mixture 07 seems to be 
a crossover, conatus structure, sharing a substantial amount of disorder with the presence of three-fold coordi- 
nated "unit cells"; the latter can be seen as vestigial in fluid system 06. Fluid systems 09 and 10 show that this 
sponge<->gyroid structural order-disorder transition not only occurs via an increase in surfactant concentration 
(a lyotropic transition), but in the interaction strength between surfactant with itself and with the interface. 
We leave for further work a systematic investigation of the {n^^^^, gbs, 5ssj ffbr} parameter space in mapping 
out the equilibrium mesostructures' phase diagram. In this endeavour, recently developed compusteering tools 
|53l Is^ may prove valuable in optimising expensive simulation time: they allow the user to postprocess and 
visualise the compute job's output at run-time with negligible turnaround times, and eventually temporarily 
stop execution in order to modify simulation parameters which are fed back into the algorithm on immediate 
restart. 
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Finite size effects can play an important role in the stabilisation of fluidic structures like these, given that we 
are using periodic boundary conditions. With this in mind, and using the same parameters as for mixture 09, 
we computed the wavenumber-averaged difference (Aj^^n'S) for each time step of evolution of the spherically 
averaged structure function S{k) between lattices of sizes and N'^, where N,N' = 64, 128,256. Note that 
the lattice size is increased eight and sixty-four times from the original 64"^ size. Finite size effects would be 
present if (Aatjv'S') were larger than the error derived from the differences and the averages. Nonetheless we 
found (An^n'S) to be larger than the error (27% larger on average for N = 128 and N' = 256), the fact that 
it strongly decreased with N (i.e. (Ai28,256'S') ~ 0.38(A64,i28'5')) provides the confidence necessary to assert 
that finite size effects are not significant in the N = 128 simulations we are about to report. Moreover, as we 
shall see immediately, since the structures corresponding to a 64'^ lattice exhibit the same morphologies as do 
the 128'^ and 256'^ cases, the qualitative features of the former are the same as those for the asymptotic limit 
N oo. This also extends to the oscillation of the structure function: the equivalent of Fig.|Slfor the 128^ and 
256'^ cases (not shown) exhibits similar features, albeit including more wavenumbers for which S{k) grows. 

Figure El displays three viewpoints of the isosurface (j) = 0.40 for the same composition of fluid mixture 09 as 
simulated on a 128^ lattice. We show the restriction to a 33^ subdomain, all at time step t — 15000, together 
with the oil-water interface. Whereas on 64^ (or smaller) lattices the liquid crystalline structure uniformly 
pervades the simulation cell, on 128'^ (and larger) lattices there are some imperfections present resulting in 
ordered subdomains with slightly varying orientations between which there exist domain boundaries. These 
boundaries can be considered as defects in the structure, the presence of which is a characteristic feature of 
liquid crystals. A time scale for the dynamics of some of these defects for our simulated gyroids (simulation 
run 09) can be roughly estimated from Fig. E| we observe for that particular slice that the topological genus 
of the interface changes in an interval ranging between 500 and 1000 time steps. State-of-the-art visualisation 
proved key in the analysis of results, and virtual reality technologies can enhance its usefulness by increasing 
interactivity with the data |[5Jj. 

Small-angle x-ray scattering (SAXS) techniques have been widely used in the determination of the nanostruc- 
ture of fluid mesophases 13, '53, Issl Isq . SAXS spectra, or their numerically computed versions |53| , give peak 
patterns for these mesophases that are used as fingerprints in determining unknown structures. However, the 
lattice resolution of our simulations is insufficient to detect multiple peak fingerprints in plots of the spherically 
averaged structure function. Instead, its unaveraged counterpart, S'(x), shows complete agreement of ratios of 
reciprocal vector moduli with those observed in diffraction patterns of the gyroid, as we display in Fig. 1131 for 
fluid composition 09 . In addition, visual inspection of the unit cell of the oil- water interface unequivocally 
identifies it with that of the gyroid. The size of such a unit cell as seen in optical textures allows us to asso- 
ciate a length scale to the lattice for a particular experimental realisation. For the system reported by Hajduk 
et al. |55l |. the lattice would need to be 291 nm in side length with a resolution of 2.3 nm per lattice unit. 

Although previous simulation papers on amphiphilic mixtures using free-ene rgy based Langevin diffusion 
equations have reported the reproduction of structures resembling the gyroid [sj, [s^l, none of them have 
studied its features or dynamics, or incontestably demonstrated its gyroid morphology. Furthermore, in one of 
these articles we observe that the fluid mesostructure is not stationary 31], whereas by eyeballing the whole 
simulation cell in another |3^ one becomes aware that the structure has a morphology which is reminiscent of 
the molten gyroid we describe here. 

In our simulations we observe that, at late times, the gyroid is much closer to stationarity than the sponge 
mesophase: for equal time slices in their evolution, temporal changes of the mesostructure over a period of 
1000 time steps are considerable for the sponge (e.g. simulation run 06) whereas for the gyroid (e.g. simulation 
run 09) they appear as slight interfacial rearrangements and undulations, reminiscent of breathing modes, 
keeping the variation in the position of each unit cell small compared to the lattice size N . This late-time 
(ca. 30000 time steps) structural dynamics is characterised in our simulations by the fact that the topological 
genus is (statistically) preserved for the gyroid; for the sponge, it is not. This can be understood as structural 
stabilisation by rigidity in the gyroid, and a flowing, glassy dynamics for the sponge. Such a distinctive behaviour 
for the sponge may have a bearing on its density fluctuations and render them different to those occurring in 
a topology-conserving dynamics. It is therefore not surprising that (a) we found the oscillation modes for the 
sponge (simulation run 06) to be at least one order of magnitude more intense than those for the different 
gyroids we simulated (simulation runs 08, 09 and 10, c/. Fig. |3, and (b) recent experimental studies, using 
dynamic light scattering and various relaxation methods |58j . do not report on fluctuations for the gyroid j59l |. 
whereas they do for the sponge mesophase. 

It is accepted wisdom [60j | , and a working hypothesis in many simulation studies psl l6lL l62j , that periodically 
modulated phases may arise in fluid mixtures whenever a repulsive long-range interaction competes with an 
attractive short-range one for a conflguration that minimises the interfacial Hamiltonian, possibly also in the 
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presence of a thermal, entropic contribution. Little is said about whether non-locality is not only a sufficient 
but also a necessary condition, or whether the non-locality of the relevant model needs to be imposed ab initio 
or is rather an effective emergent feature picked up by the order-parameter autocorrelation function. The LB 
model we employ in this paper only incorporates local interactions in its mesodynamics; this feature allows 
its algorithm to be easily parallelised and achieve exceptionally good performance |^J. Nonetheless, we have 
demonstrated that the model is able to simulate liquid crystalline, cubic mesophases, such as the gyroid in binary 
immiscible fluid mixtures with an amphiphilc, and the primitive "P" in binary, amphiphile-solvent mixtures 
|24j, whose density-density correlations are markedly non-local, and in the formation of which hydrodynamic 
interactions play a vital role. Non-locality, in our case, is an emergent property of a local model. 

It is worthwhile pointing out that Prinsen et al. |23| , using a DPD model and basing their claims on Monte 
Carlo studies of equilibrium cubic phases by Larson |63|, suggested that cubic phases could be engineered 
to appear if their bead-rod model were elaborated beyond dimers. If fact, Groot & Madden's (inconclusive) 
finding of a gyroid-like structure with a DPD model for bead-spring chain copolymers |2l| might be considered 
to support such an assertion. Part of the importance of our simulations, as well as those of Nekovee & Coveney 
|24| . is to refute such conjectures, by demonstrating that cubic mesophases arise in very simple, minimalist, 
athermal and hydrodynamically-correct fluid models, with a locally-interacting vector order parameter and 
reproducing universal behaviour. 

VII. SUMMARY AND CONCLUSIONS 

Our simulations furnish the first quantitative account of phase segregation kinetics and mesophase self- 
assembly in amphiphilic fluids with a three-dimensional model based on the Boltzmann transport equation. 
The method is hydrodynamically correct, athermal, and models the amphiphilic species as bipolar, point-like 
particles experiencing short-range interactions with mean fields created by the surrounding binary immiscible 
( "oil- water" ) and amphiphilic medium. 

We studied the phase segregation pathway in a homogeneous oil-water-surfactant mixture at composition 
1 : 1 : X, respectively, where < a; < 1.3 is the surfactant-to- water (or to oil) mass fraction. We observed 
segregation slowdown in the average size of oil-water domains with increasing x, and the reproduction of known 
crossovers, namely, from algebraic to logarithmic to stretched exponential functions. This confirms the usefulness 
of our method in apprehending the fundamental phenomenology of amphiphilic fluid mixtures; the presence of 
transients in these crossovers is gratifying given their experimental observation. In order to rule out an increase 
in total density as x is increased as a factor contributing to the slowdown along with the reduction in surface 
tension, future work should investigate domain growth at constant total density. 

The stretched exponential functional form occurring at domain growth arrested regimes can be ascribed to 
the accumulation of a large number of relaxation modes associated with surfactant dynamics onto and at the 
oil- water interfaces [T9j |. The late-time structure at these regimes are (a) disordered non-stationary sponge 
mesophases if the initial amphiphile concentration is lower than a threshold region or (b) well-defined liquid- 
crystalline cubic gyroid mesophases of pinned domain sizes with defects if this initial concentration is higher 
than such a threshold. We also found thresholds in the surfactant coupling strengths, \gss\ and \gbs\, for a 
sponge-to-gyroid transition. In the transition region we observed a crossover structure sharing the structural 
features of both the gyroid and sponge. We also found that, for the number of time steps simulated, both sponge 
and gyroid exhibit undamped oscillations at all length scales. For some length scales, the temporal trend of 
their Fourier amplitudes is to slowly die out; for others, it increases. 

We found that extremely slow domain growth can be mistaken for genuine arrested growth if attention 
is not paid to minute length scales. Truly segregation-halted regimes exhibit oscillations in average domain 
size, which can be seen at sufficiently late times. These oscillations are caused partly by Marangoni flows 
generated by inhomogeneities in the surfactant adsorbed on the oil-water interface, and partly by a surfactant 
dynamics dictated by competing mechanisms, namely, surfactant attraction towards the interface and surfactant- 
surfactant interactions. Because our model does not presuppose that all the surfactant is adsorbed on the 
interface, as Langevin approaches based on the adiabatic approximation do psl l28l l60l| . surfactant-surfactant 
interactions are not limited to repulsion. Hence, in regimes of large surfactant concentration, and especially 
in those for which regions of the (oil-water) interface can be sufficiently close to each other, surfactant is not 
constrained to dwell on the interface; rather, it is reasonable to propose the existence of an adsorption-desorption 
dynamics driving surfactant towards and away from it. Our results showing that (a) an increase in \gss\ excites 
higher frequencies, (b) an increase in \ghs\ dampens most frequencies, and (c) there appear surfactant currents 
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bridging adjacent interfacial regions, confirm this proposal. 

Our method is not only the first lattice-Boltzmann model to deal with segregation kinetics in three-dimensional 
amphiphilic fluid mixtures, but the first complex fluid model to unequivocally reproduce the gyroid cubic 
mesophase, using a high level of abstraction in modelling the amphiphile. The truly mesoscopic, particulate 
nature of the surfactant in this model accounts for the complex, dynamical behaviour observed, even in a noiseless 
scenario like ours. It is not surprising that Ginzburg-Landau-based Langevin models which treat surfactant 
implicitly through a scalar parameter modifying the free energy are only able to exhibit oscillations which decay 
rapidly in time and whose frequency spectrum has but one peak. Our simulations of liquid crystalline mesophases 
prove that, contrary to what has been previously claimed |63,|63,|63|, surfactant-surfactant interactions need not 
be long ranged in order for periodi cally modulated, long-range ordered structures to self- assemble. In addition, 
our findings rebut the suggestion |2a | that cubic mesophases can be simulated only when the amphiphile is 
modelled with a high degree of molecular specificity. 

The simulation of the sponge and gyroid phases, and their complex oscillatory dynamics, confirms the richness 
of our model's parameter space. Our lattice-Boltzmann model provides a kinetically and hydrodynamically 
correct, bottom-up, mesoscale description of the generic behaviour of amphiphilic fluids, which is also extremely 
computationally efhcient on massively parallel platforms. Future extensions of this work include the search for 
regimes leading to equilibrium mesophases of more varied symmetries, the study of shear-induced symmetry 
transitions, and large scale studies of defect dynamics in liquid crystalline phases. In fact, the TeraGyroid 
project, a successful Grid-based transatlantic endeavour employing more than 6000 processors and 17 teraflops 
at six supercomputing facilities 0|, has its scientific raison d'etre based on the results we report in this paper. 
TeraGyroid proves the value of computational steering tools '53^ in mapping new parameter space regions of 
our model. 
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FIG. 1: Surface tension dependence on the surfactant concentration (mass fraction, cf. Tabled as measured at a planar 
interface making use of Eq. 1)15^ . A lattice of size 4 x 4 x 128 was allowed to evolve up to time step 25000, and pressure 
tensor components were measured every 1000 time steps. The surface tension tends to grow with time and reaches a 
horizontal asymptote; at that time step the surface tension only differs in 16% from that at the previous measurement. 
Interpolation serves as a reference to the eye. Coupling constants used were pbr ~ 0.08, ^bs = —0.006, and g^s = —0.003. 
Oil and water densities used were n'"-"^ = n'-'^'^ = 0.7. All quantities are reported in lattice units. 
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FIG. 2: Temporal evolution of the average fluid-fluid domain size for surfactant concentrations 0.0, 0.15, 0.22, 0.30, 
0.35, 0.40, 0.60 and 0.90 for curves from top to bottom and corresponding to simulation runs 01, 02, 03, 04, 05, 06, 
07 and 08, respectively (c/. Table HJ. Measurements have been taken every 25 time steps, and the plots include error 
bars, which represent the uncertainty (one standard error) transmitted from the standard error of the structure function 
spherical average. We used a lattice of size 128'^ for simulation run 01 and 64^ for the remaining curves, since finite size 
effects start to creep in for domain sizes larger than L ~ 30. All quantities are reported in lattice units. Note that the 
surfactant-containing fluids lack the zero-growth, linear transient found for simulation run 01 flU . 



19 




50 r— 
40 - 

m- 




10 ' 

I \ \ \ 

7 7.5 S 8.5 

In I 

(b) 

FIG. 3: Panel (a) shows the time evolution of the average domain size for simulation runs 01, 02 and 03, see Fig. |5| 
The log-log scale helps to visually detect behaviours following Eq. 11811 — in this case, that of the uppermost curve. The 
straight line above it serves as a guide to the eye only and its slope is given by ci in Table |n] Panel (b) shows the 
evolution of the average domain size with the logarithm of the time step for simulation runs 01, 02 and 03, on a log-log 
scale. This is useful to discriminate growth between that of Eq. p8|l and Eq. I|19|l : see Table ITT| for the fitting parameters. 
The straight solid line shown indicates a good fit to Eq. 1191 for t > 1100 {\nt ~ 7). For t < 1100, the curve is better 
fit by Eq. I|18^ . albeit still quite poorly. Measurements have been taken every 25 time steps, and the plot includes error 
bars representing the uncertainty (one standard deviation) of the spherical averaged structure function. We use a 128^ 
lattice for simulation run 01 and a 64'^ lattice for the remainder. All quantities are reported in lattice units. 
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FIG. 4: Spherically averaged structure functions for the oil-water order parameter simulation run 06 (c/. Table UTTl . 
According to how close to asymptotic behaviour the distribution of domain sizes appears to be, we have classified 
simulation times for this simulation run in three groups: early times (time steps 25, 50, 75, 125, 150, 175, and 300 in 
the plot), intermediate times (time steps roughly from 800 to 1700), and late times (time steps 1800, 2300, 2800, 3300, 
3800, 4300, 4800, 6000, 10000, 14000, 18000, 22000, 26000, and 30000 in the plot). Error bars represent the standard 
error of the shell average. Lattice size is 64"^. All quantities are reported in lattice units. 
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FIG. 5: Log-linear plots of the spherically averaged structure functions at time step 7200 for increasing surfactant 
concentrations indicated by the numerical labelling on each curve, corresponding to simulation runs 01, 02, 03, 04, 05, 
06, 07 and 08, respectively, c/. Table Q] For large wavelengths, panel (a), we can see how the peaks move to higher 
wavenumbers, decrease in height, and broaden. Note that for short wavelengths, panel (b), the only straight tail is for 
curve n'"''^ — 0.0, whose slope is —4.46 x 10~*. Error bars represent one standard error of the shell average S{k). Lattice 
size is 128'' for simulation run 01 and 64^ for the others. All quantities are reported in lattice units. 
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FIG. 6: Temporal evolution of the average domain size for simulation runs 06, 07, 09, 10, and 08, as seen from top 
to bottom at t = 10000 (c/. Table ITTTt . Measurements have been taken every 25 time steps; error bars are included 
and represent the uncertainty ( "one sigma" ) transmitted from the standard error of the spherically averaged structure 
function. Caveat lector, an oscillation in the average domain size is genuinely representative of oscillations in the domain 
sizes only if error bars are smaller than the oscillation amplitude. Lattice size is 64"^. All quantities are reported in lattice 
units. 
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FIG. 8: Temporal dependence of the structure function for simulation run 09, c/. Table Uni Panels (a) and (b) show 
the short and long wavelengths, respectively. Measurements have been taken every 25 time steps; error bars have been 
included. Lattice size is 64"^. All quantities are reported in lattice units. 
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FIG. 9: Frequency power spectra of the structure function for simulation runs 06, 07, 09 and 10 (c/. Table HID at 
wavenumber 0.589. Panel (b) is a blowup of panel (a) for long periods of oscillation. Error bars were neither considered 
for the Fourier analysis nor plotted here. All quantities are reported in lattice units, and uj is the inverse period. 
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FIG. 10: Slices ol the surfactant density in 33'^ subsets of a 128^ lattice, for composition 09 (c/. Table^l. Panels (a), (b) 
and (c) are snapshots at t = 14000, 14600 and 15000 time steps, respectively, which are times for which the structure 
is close to equilibration. We only show regions where surfactant density is the highest (0.31 < < 0.35), in grey and 
white, where lighter shading denotes a higher value. We can see that the surfactant mainly concentrates around the 
oil-water interface ((?i(x) = 0), whose intersection with the slice is depicted as open undulating or closed lines. Also, 
there are ordered, crystalline regions along with smaller regions lacking long-range order and evolving in time. Finally, 
is non-uniform on the interface (as regions with lighter shading show) , favouring the formation of "surfactant bridges" 
between adjacent portions of the interface; this leads to Marangoni effects which account for the observed oscillatory 
behaviour. All quantities are in lattice units. 
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FIG. 11: Equilibrium minority- phase order parameter isosurfaces (^(x) = 0.37 (in lattice units), taken at time step 
t = 30000, at which —0.69 < 4> < 0.68 over the whole lattice. Snapshot (a) corresponds to simulation run 06, snapshot 
(b) to simulation run 07, and snapshot (c) to simulation run 08, c/. Tabled Shown are 16 < 2 < 48 slabs of a 64 
lattice. Note that increasing surfactant concentration leads to an increased ordering in the mesostructure: simulation 
run 06 exhibits sponge-like features, whereas simulation run 08 is a liquid crystalline cubic gyroid mesophase. Snapshot 
(b) is a crossover state in this lyotropic transition — a "molten gyroid" . 
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FIG. 12: Isosurfaccs of the order parameter at time stop t = 15000 in a highly ordered 33^ subdoinain on a 128^ lattice for 
fluid composition 09. Panels (a), (b) and (c) display the (/'(x) = 0.40 (in lattice units) minority phase isosurface viewed 
as axononietric projections along the (10 0), (1 1 1) and (1 10) directions, respectively. Panel (d) shows the interface of 
the same lattice subdomain along direction (111), where black and white have been used to distinguish one immiscible 
fluid phase from the other. 
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FIG. 13: Temporal evolution of {kx,ky) slices of the structure function S(k) for fluid composition 09, viewed along the 
(10 0) direction. Panels (a) and (b) show slices at time step t = 500, where kz/{2n/N) — for (a) and kz/ {2tt/N) = ±14 
for (b). Similarly, panels (c) and (d) are slices at time step t = 15000 for kz/{2n/N) = 0, ±14, respectively; A'' = 128 
is the lateral lattice size. Shading denote intensities S = 1,50 and 100, where lighter greys up to white mean higher 
intensities. The spherical shell structure in (a) and (b) indicates the presence of a sponge (microemulsion) phase, which 
becomes anisotropic at later times, (c) and (d). The superposition of slices (c) and (d), namely, the ratio of peaks' 
positions of S'(k), are in full agreement with SAXS experimental data for the gyroid mesophase. All quantities are 
in lattice units. 
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